Density fluctuations and the structure of a nonuniform hard sphere fluid 
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We derive an exact equation for density changes induced by a general external field that corrects 
the hydrostatic approximation where the local value of the field is adsorbed into a modified chemical 
potential. Using linear response theory to relate density changes self-consistently in different regions 
of space, we arrive at an integral equation for a hard sphere fiuid that is exact in the limit of a slowly 
varying field or at low density and reduces to the accurate Percus-Yevick equation for a hard core 
field. This and related equations give accurate results for a wide variety of fields. 
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Determining the structure and thermodynamics of a 
hard sphere fluid in a general external field is a basic prob- 
lem in the theory of nonuniform fluids . When the field 
represents walls, slits, or pores, this is the simplest model 
system describing the interplay between excluded volume 
correlations in the bulk fluid and the effects of confining 
geometries. Other choices for the field can represent fixed 
fluid particles and thus give information about multiparti- 
cle correlation functions (2|J^. Moreover, within mean field 
theory, a hard sphere fluid in the presence of an appropri- 
ately chosen "molecular field" potential can also describe a 
liquid- vapor interface and wetting and drying phenomena 
near walls or other solutes for simple liquids 

There are two limits where accurate solutions to this 
problem are already known. The simplest limit occurs 
when the external field (f){r) varies so slowly that it is es- 
sentially constant over the range of a correlation length in 
the bulk hard sphere fluid 0| . The partition function and 
the density p{r; [(/>], /i^) in the grand ensemble are func- 
tional of the external field and functions of the chemical 
potential /z^, and depend only on the difference between 
these quantities |^]. Thus we can subtract any constant 
from both /i^ and (p with no effect on the exact structural 
or thermodynamic properties of the fluid. In particular, 
for any fixed position ri we can define a shifted field 
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whose parametric dependence on ri is denoted by a 
superscript, and we have for all r the exact relation 
p(r; [0], /x^) = p(r; [(/('"ijj/x'"!). However, by construction 
the shifted field 0""^ (r) vanishes at r = ri and it remains 
very small for r near ri when is very slowly varying. In 
such a case we have 
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Here p{fJ-^^) = p^^ is the hydrostatic density: the density 
of the uniform fluid at the shifted chemical potential /i''^ . 



From Eq. (g) this depends only on the local value of the 
field (/) at ri . Equation (||) is exact in the limit of a very 
slowly varying </>, even when its magnitude is very large 0] . 
However Eq. (||) must fail if cj) has significant gradients over 
the range of a correlation length in the bulk fluid. 

Accurate results can also be found in the limit of a hard 
core field, where varies as rapidly as possible. Then 
the Percus-Yevick (PY) equation is often very accurate, 
particularly at low to moderate densities Simple cor- 
rections to the PY approximation for hard core fields such 
as the generalized mean spherical approximation (GMSA) 
are available |^ and give even more accurate results. 

Recent work |9|J^ has provided a new interpretation of 
the PY approximation for hard sphere fiuids that is physi- 
cally suggestive and suitable for generalization. Computer 
simulations have shown that even large spontaneous den- 
sity fluctuations in a uniform hard sphere fiuid can be 
accurately described using the same Gaussian probability 
distribution that controls small fluctuations pO[ |. Small 
density changes induced by small changes in the fleld 
for a general system with chemical potential p,^ , temper- 
ature ksT = f3~^, and density p{r;[(j)], p^) = p{r) are 
linearly related: 
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through the linear response function 

X"^(ri,r2; [p]) = J(ri-r2)/p(ri)-c(ri,r2; [p]). (5) 

Here c is the direct correlation function of the system; it is 
a functional of the density p{r). Note that the external po- 
tential appears explicitly only on the left side of Eq. (^). In 
most standard applications, one considers perturbations 
about = 0, so p{r) — p^ and reduces to the uniform 
fluid function x~^(''i2; P^)- By the fluctuation-dissipation 
theorem, this same function controls the small Gaussian 
density fiuctuations in the uniform fiuid |^,^. 

Can Eq. also be used to calculate the response to a 
large (hard core) field? Certainly this linear relation be- 
tween a finite external field perturbation and the induced 
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density must fail for values of ri where the field is very 
large. Conversely, Eq. (||) should be most accurate for 
those values of ri where the field is small — in particular 
where the field vanishes, and the simulations suggest that 
even large fluctuations in the absence of a field are well 
described by a linear fluctuation theory. 

If we use Eq. (Q) only for those values of ri where the 
field perturbation vanishes, then through the integration 
over all r2 it relates density changes in the region where 
the field vanishes to the density change in the region where 
the field is nonzero |6|. For a hard core perturbation the 
density in the latter (hard core) region must vanish exactly 
and we can use Eq. (jj) to determine the induced density 
changes in the region of zero field outside the core if the 
X~^(ri2; p^) for the uniform fluid is known. 

For a general hard core solute (e.g., a hard wall), this is 
equivalent to the FY approximation for the solute-particle 
direct correlation function An equation equivalent 

to the FY equation for hard spheres results from a self- 
consistent application of this procedure for a potential rep- 
resenting a hard sphere fixed at the origin. By making use 
of the exact relation between the density induced by such 
a fixed particle and the pair correlation function of the 
uniform hard sphere fluid |^ , one obtains the FY approx- 
imations for both the induced density and X~^{^12', P^)- 

Frogress was possible in this special case because we 
could impose the exact "core condition" that the 

density vanishes in the hard core region. But what can 
be done for a general finite external field where the asso- 
ciated density response is not known in advance? 

We note that there is a common feature of both lim- 
its discussed above. Simple and accurate approximations 
are available wherever the local value of the field vanishes. 
Thus depending on the value of the external field we can 
apply the exact hydrostatic shift in Eqs. (|l|) and (||) at 
each point of space to ensure that this optimal condition 
holds true. The hydrostatic limit for a slowly varying field 
is automatically satisfied. At the same time we can use 
accurate methods based on linear response theory that re- 
late only densities in different parts of space to take into 
account nonlocal effects of the external field. In contrast, 
standard methods generally use only the response function 
of the bulk fluid and try to prescribe nonlinear closures 
that directly relate the field and the density [||. 

To develop a quantitative theory we first note that if 
(j) depends on a parameter A (denoted by the asso- 
ciated density is p\) then Eq. (^ can be rewritten ex- 
actly by setting S(f>{ri) — [d(j)x{ri)/dX]dX and Sp{r2) — 
[dpx(r2) / dX]dX. In view of the discussion above we want 
to remove the explicit appearance of the field on the left 
side of Eq. (^ . To that end we rewrite Eq. (^) using the 
shifted field and chemical potential of Eqs. (^ and (|^). We 
then introduce a modified potential (//^ (r) depending on a 
coupling parameter A with < A < 1 such that d(tf.^ (r) / dX 
vanishes for all A at r = ri; for A = 0, vanishes every- 
where and the associated density p'^^ (r) = p'^^ , the uniform 



hydrostatic density at ri, while for A = 1, (jf^ = (jf'^ and 
the associated density is p(r; [cff^],^^^) = p{r), the desired 
density of the fluid in the general external field 0. One pos- 
sible choice is (j/^^ (r) = A^""^ (r), but the linear dependence 
on A is not essential in the following. 

Thus we find = Jdr2X'Hri,r2;[pl'])dp\'{r2)/dX 
from Eq. (|^) and we can integrate with respect to A to 
obtain the formally exact result: 

= Jdr2 dX x"'(ri,r2; [p\']) dp^^ {r2) / dX . (6) 

This equation has the desired features that the external 
field does not appear explicitly and the induced density 
changes in different parts of space are related to one an- 
other through a linear response function. 

For practical calculations we must approximately carry 
out the A integration. The simplest treatment approxi- 
mates x~^(ri, [p]^]) for all A by X~^(''i2; its value 
at A = 0, consistent with the idea that even large changes 
in the density are controlled by the same response func- 
tion as in the case of small changes. The A integration 
involving the density can then be carried out exactly, and 
using Eq. (^ we obtain our final result, the hydrostatic 
linear response (HLR) equation: 

p(ri) = p^- + p^- j dr2c(ri2;p'-^)[p(r2) - p^']. (7) 

Equation (^ is a linear integral equation relating the 
density p{yi) at a given ri to an integral involving the 
density p{y2) at all other points and a uniform fluid ker- 
nel c{ri2\p^^) that depends implicitly on ri through p^^ . 
This new feature presents no technical difficulties in de- 
termining a self-consistent numerical solution. We found 
that Ficard iteration works very well. See Refs. (TTj and 
for details about the numerical solution. 

Equation (Q) has the following remarkable properties, i) 
It is exact when 0(r) is very slowly varying, ii) It is exact 
for any (j){v) at low enough density, where there is a local 
relation between the potential and induced density, iii) 
For a field (/)(r) from a general hard core solute, Eq. (Q) 
reduces to the FY approximation, as discussed above. Any 
desired representation of the uniform fluid c can be used. 

Other equations can be derived by making different 
approximations while carrying out the A integration in 
Eq. (||). We note from Eq. (^) that the local (5— function) 
part of becomes relatively more important as the den- 
sity yo(ri) tends to zero either because of a harshly repul- 
sive (/)(ri) or because p^ is small. A somewhat inconsistent 
approximation that does however exactly describe the lo- 
cal part of x"^ sets x"^(ri,r2; [p^^]) ~ (5(ri-r2)/yO^' (ri)- 
c{.fi2] p'^^), thus keeping the exact A-dependence in the lo- 
cal part of x~^ but setting A = in the nonlocal part. 
The A integration can then be carried out exactly and we 
arrive at an alternative hydrostatic mixed (HM) equation: 
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FIG. 1. Correlation functions for hard spheres in the pres- 
ence of spherical parabolic potentials shown in the inset (solid 
lines) as given by theory and simulation. The upper curve 
corresponds to the smaller potential and has been displaced 
upward by one unit. Also shown in the inset (crosses) are the 
potentials predicted by Eq. (M) given the simulation data. 
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FIG. 2. Correlation functions arising from the steep planar 
triangular potential shown in the inset for = 0.55 (lower 
curve) and = 0.75; the latter curve has been displaced 
upwards by one unit. 



p(ri) = p^- exp{ j dr2c(ri2;/0[p(r2) - p"']}- (8) 

In regions where the external field 0(ri) is infinite both 
the hydrostatic density p""^ and the exact density p{yi) 
must vanish; the solutions to both Eqs. (0) and (||) clearly 
satisfy this condition. Moreover, the p(ri) solving Eq. (ph 
is always nonnegative, which is not the case for Eq. (M), 
and properties i) and ii) above still hold true. As we will 
see, Eq. (|^) also gives very good results for a wide variety 
of potentials. However, for a hard sphere potential Eq. (||) 
reduces to the hypernetted chain (HNC) equation for 
hard spheres, which is known to be much less accurate 
than the PY equation at high density 

To give some indications of the accuracy of Eqs. (0) 
and (|g) we report solutions for a number of different ex- 
ternal fields and compare with the results of computer 
simulations [^. We first consider two model potentials 
designed to show both the strengths and the weaknesses 
of the present methods. Then we consider more realistic 
potentials arising from a mean field treatment of wetting 
and drying phenomena in the Lennard Jones (LJ) fluid. 

Fig. 1 shows the correlation function g{r) = p(r)/ for 
a hard sphere system at a moderately high bulk density 
p^ = 0.49 in the presence of two deep attractive spher- 
ical parabolic model potentials shown in the inset. (Re- 
duced units, with distances measured in units of the hard 
sphere diameter are used.) Both equations reproduce the 
increased density inside the well, and the nonlocal oscil- 
latory excluded volume correlations, which show a local 
density minimum at the center of the well due to packing 
effects. 

Since the external field enters Eqs. and (|^) only lo- 
cally through its effect on p""^ , it is also easy to use these 
equations for the inverse problem of determining the field 



associated with a given density profile. As an example, 
the crosses in the inset gives the potentials predicted by 
Eq. (^ given the simulation data for g{r). 

One would expect methods based on expanding about 
the uniform hydrostatic fluid to be least accurate for po- 
tentials with very steep gradients, exemplified by the re- 
pulsive planar triangular barrier potential shown in the 
inset to Fig. 2. Even here reasonably good results are 
found from Eqs. (|^) and (^) at the moderately high den- 
sity of p^ — 0.55 (with excellent results at lower densities). 
However, at p^ = 0.75 very noticeable errors are seen in 
the shape and height of the first peak and the amplitude 
and phase of subsequent peaks in both approximations. 

Results actually improve as the barrier height increases 
and the potential approaches a hard wall potential: the 
theory does better for hard cores because there is no re- 
gion of space where there is a large gradient in the external 
field while at the same time the local density is nonzero. 
Equation (|^) then reduces to the accurate hard core PY 
wall-particle equation, and Eq. (||) describes subsequent 
peaks better, though the characteristic HNC overshoot of 
the height of the first peak for hard core systems becomes 
very evident as the slope increases. Large repulsive poten- 
tials with very steep gradients are better treated by "blip 
function" methods Hj or other expansions about a hard 
core system. 

In most realistic applications we can envision, the ex- 
ternal potential can be divided into a harshly repulsive, 
essentially hard core, region well treated by the present 
methods (with "blip function" corrections for finite core 
softness), and an extended interaction region where the 
potential varies slowly enough that the theory again gives 
good results. This is the basis of the two step method used 
to solve the mean field equations for the LJ fluid in Refs. 
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FIG. 3. Correlation functions (vertically displaced by 0.1 
unit) induced by the potentials shown in the insets. These 
potentials arise from a generalized mean field treatment of the 
response of a LJ fluid to a hard sphere solute with varying 
size. Circles indicate simulation data and the solid line give 
the predictions of Eq. . 



for a variety of different problems. 

In Fig. 3 we show the reference system g{r) as predicted 
by the HLR equation induced by the repulsive external 
fields shown in the inset. These fields arise from a gen- 
eralized mean field treatment describing the density 
distribution in a L J fluid for a state very near liquid- vapor 
coexistence with a reduced density of 0.70 and reduced 
temperature of 0.85 as the radius of a hard sphere solute 
is varied . Units are reduced by the usual L J potential 
parameters a and e. The upper curve shows a molecu- 
larly sized solute, which induces density oscillations like 
those seen in the radial distribution function. For larger 
excluded volume regions, the unbalanced attractive forces 
Q in the LJ fluid contribute a stronger and longer ranged 
repulsive component to the effective potential in the as- 
sociated hard sphere system as seen in the inset and par- 
tial drying occurs. This manifests itself in the middle and 
lower curves by smoother profiles, with a density minimum 
near the solute. Other limits of this problem, such as hard 
spheres near a hard wall or the mean field approximation 
to the smooth vapor-liquid interface of a LJ fluid, are well 
described using the same HLR equation in the presence of 
the appropriate effective mean field. 

While we have concentrated on the structure of 
hard sphere systems here, even large density fluctuations 
in more complicated liquids like water are also well de- 
scribed by a Gaussian model in the presence of 
an effective field [|l9| describing the unbalanced attractive 
forces arising when interfaces form. By taking this field 
into account, Lum, Chandler and Weeks jl^] were able 
to extend the (field free) Gaussian theory of hydropho- 
bicity of Pratt and Chandler for small hydrophobic 
solutes to large hydrophobic solutes, where drying inter- 



faces, very similar to those seen in Fig. 3, are predicted to 
occur. More generally, we believe that extensions of linear 
response methods or Gaussian fluctuation theory |2^] to 
include the effects of appropriately chosen external fields, 
with optimal treatment of the field by equations similar to 
the HLR equation (^, will prove useful in a variety of dif- 
ferent problems. This work was supported by the National 
Science Foundation, Grant No. CHE9528915. 
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